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A Jastrow wave function (JWF) and a shadow wave function (SWF) describe 
a quantum solid with Bose-Einstein condensate; i.e. a supersolid. It is 
known that both JWF and SWF describe a quantum solid with also a finite 
equilibrium concentration of vacancies x v . We outline a route for estimating 
x v by exploiting the existing formal equivalence between the absolute square 
of the ground state wave function and the Boltzmann weight of a classical 
solid. We compute x v for the quantum solids described by JWF and SWF 
employing very accurate numerical techniques. For JWF we find a very small 
value for the zero point vacancy concentration, x v = (1.4 ± 0.1) x 10~ 6 . For 
SWF, which presently gives the best variational description of solid 4 He, we 
find the significantly larger value x v = (1.4 ± 0.1) x 10~ 3 at a density close 
to melting. We also study two and three vacancies. We find that there is a 
strong short range attraction but the vacancies do not form a bound state. 



1. INTRODUCTION 

A fascinating topic in modern quantum many-body Physics is the pos- 
sibility of a supersolid state of matter, i.e. a solid which displays some 
superfluid properties. Since the first theoretical speculations on its exis- 
tence in the early seventies, supersolid has gained a wide attention.^ Being 
such a supersolid phase a direct consequence of macroscopic manifestation 
of quantum properties, the highly quantum solid 4 He was soon recognized 
to be the natural candidate to realize it; but it has revealed to be elusive to 
the experimental observation until 2004, when the observation of non clas- 
sical rotational inertia effects (NCRI) 2 provided the first possible signature 
of its presence. Many experiments followed the original onepEl giving rise 
to a puzzling scenario: NCRI is now observed in torsional oscillator exper- 
iments both when solid 4 He is inside a porous material as well as in bulk 
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solid He, but the detected supersolid fraction p s shows a strong dependence 
on the details of the experimental realization. The observed dependence of 
p s on the sample, on 3 He concentration^ and also on annealing } 2 ^ 5 ! leads 
to a largely diffuse conviction that NCRI should be essentially ascribed to 
extrinsic disorder rather than to an intrinsic property of solid 4 He. Some 
aspects of the measurements (as the metastability at low temperatures, the 
value of the critical velocity, p s which vanishes in a continuous way with a 
dissipation peak, the dependence of p s on the oscillator frequency) seem to 
point in the direction of a liquid of vortices as suggested by AndersonP On 
the other hand, recent experiments with an high quality 4 He crystal and 
the observation of a peak in the specific heat in correspondence of the super- 
solid transition® seem to put some restrains to the interpretation of NCRI 
as purely extrinsic effect. There is no doubt that defects play a relevant role 
in NCRI, but what happens if a more and more perfect crystal is grown? 
Will the supersolid response disappear or some NCRI effect survives? The 
answers to these questions lie in the properties of the ground state wave 
function of a quantum crystal. To what extent particles localization induced 
by the crystallization process succeeds in overcoming the zero-point motion 
of the atoms? Is it possible that even in the true equilibrium state a finite 
amount of intrinsic defects, such as vacancies {21121 may survive leading to a 
superfluid response? This ques tion c oncerning the true nature of the ground 
state of solid 4 He is still openPSHH Zero-point vacancies, i.e. vacancies in 
the ground state of solid 4 He, were the first proposed mechanism for the su- 
persolid stateP^nH Delocalized defects, acting as a dilute Bose gas, undergo 
a Bose-Einstein condensation (BEC) resulting in superfluid-like properties 
of the crystal. But zero-point vacancies have never been experimentally ob- 
served, however actual experiments cannot exclude them at a concentration 
lower than about 0.4%p3 Our understanding of the solid phase of 4 He is far 
from being complete and a full microscopic description that accounts for the 
whole supersolid phenomenology is still lacking. 

Even if it is possible to obtain essentially exact results for most prop- 
erties of solid 4 He by means of quantum simulations, e.g. PIMO 1 ^ at finite 
temperature and SPIGSp^ at zero temperature, they do not offer a simple 
straight way to answer the question of the true nature of the ground state. 
Such computations give compelling evidence that the commensurate crystal, 
i.e. a crystal in which the number of particles is equal to the number of lat- 
tice sites, has no BEC. Much more delicate is the question of the presence of 
defects like vacancies in the ground state. Whether solid 4 He, or any other 
possible quantum crystal, should have a commensurate or incommensurate 
ground state is a fundamental question due to the important role of such 
defects on the occurrence of BEC. Quantum simulation of solid 4 He with 
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one vacancy^ 17 * 18 ! 11 ! indicates that the energy of this state is above the one of 
a perfect crystal and this has been taken as evidence that no ground state 
vacancies are present in solid 4 He. However it has been argued that such 
quantum simulations present in literature do not allow to directly answer the 
question.^ In fact, in Monte Carlo simulations of crystals, the small number 
of particles, the use of periodic boundary conditions and commensuration 
effect of the simulation box with the lattice structure impose a constraint 
that makes impossible, in practice, for the system to develop an equilib- 
rium concentration of vacancies.^ This holds not only for canonical and 
micro-canonical computations (where it should be somehow obvious since 
the number of particles and the volume of the system are fixed quantities 
during the simulation), but also for grand canonical simulations.^ Further- 
more the equilibrium concentration could be so small to easily have escaped 
detection in simulations of the ground state of solid 4 HeP Therefore the 
equilibrium concentration of zero point vacancies x v can be obtained only 
indirectly, by a statistical thermodynamical analysis of an extended system, 
exactly as for classical solids 

On the other hand we know two model system that have zero point 
defects and BEC. This is the case of the quantum solid described either by 
a Jastrow wave function (JWF) or by a shadow wave function (SWF). Both 
these wave functions have a more or less successful long career as trial wave 
functions in variational descriptions of solid as well as liquid 4 He. It is well 
known that at large densities such wave functions describe a crystalline solid. 
Moreover, the presence of BEC in such solid phase is stated by a theorem 
both for JWF2H and for SWF, 1 ™ 1 and in the latter case the condensate 
fraction has recently been estimated.^ 

Furthermore, by following the Chester's argument for supersolidityf^ it 
is possible to infer that such translationally invariant wave functions (JWF 
and SWF) give rise to a finite concentration of vacancies. Even though 
this argument is common knowledge, no estimate exists of how large is such 
equilibrium concentration of vacancies, except for an approximate calcula- 
tion done some years ago for JWFP In this paper we have developed the 
route outlined by Chester and we have exploited the most advanced numer- 
ical techniques to estimate the equilibrium concentration of vacancies x v in 
the ground state of the quantum solids described by JWF and SWF. By 
a full thermodynamical analysis of the equivalent classical crystal, we have 
estimated x v = (1.4 ±0.1) x 1CT 6 for a JWF and x v = (1.4 ±0.1) x 10" 3 for 
the SWF at a density close to the melting one. 

The paper is organized as follows: in Section[2T]we introduce the Jastrow 
and the shadow wave functions briefly reporting also the Chester's argument 
for vacancy induced supersolidity; a route for the estimate of the equilibrium 
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concentration of vacancies x v is outlined in Section [23 which deals also with 
the numerical tools involved in the computation. Our result for JWF and 
SWF reported in Section |4] and discussed in Section 00 which contains 

also our conclusions. 



2. VARIATIONAL WAVE FUNCTIONS 

The ground state wave function for a system of N bosons must be 
symmetric under permutations of the coordinates, and it can be chosen as a 
real non negative function. The simplest choice for the wave function in the 
interacting bosons case is represented by the Jastrow wave function 

^j(R) = -}=\\e-^^ , (1) 

where u(r) accounts for the direct correlations among two particles and Qn 
is the normalization constant. This wave function is translationally invariant 
and describes a liquid or a solid depending on the density. The solid phase 
comes as a spontaneously broken symmetry effect at large density. Moreover 
such a wave function has been proved to have a finite Bose-Einstein conden- 
sate fraction even in the solid phasepH There is a formal identity, recognized 
long agoJ^Sl between and the normalized probability distribution in 

configurational space of N classical particles at an inverse temperature (5 
interacting via a potential 

U 3 {R) = \Y<< r v) ■ ( 2 ) 

In addition, if as commonly assumed u(r) vanishes at large distance, the 
normalization constant Qn = J dR\\ i< ^ e~ u<yTi ^ is equal to the canonical 
configurational partition function of this classical system, so that the loga- 
rithm of Qn turns out to be proportional to the excess Helmholtz free energy. 
The Chester's argument in favor of the vacancy induced supersolidity lies on 
this formal analogy.^ At large density this equivalent classical system corre- 
sponding to \ipj\ 2 is a crystalline solid. For this equivalent classical system 
the equilibrium concentration of vacancies x v (= (M — N) /M where M is the 
number of lattice sites) is non-zero even if a single vacancy has a finite cost 
of local free energy and this raises from the gain in configurational entropyP^ 
It is then straightforward to infer that also the quantum solid described by 
([1]) should have the same x v . 

One comment is in order. The pseudopotential u(r) should be deter- 
mined by minimization of the expectation value of the Hamiltonian. In the 
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liquid phase the quantitative results given by the JWF are reasonable even 
if not very accurate.^ Much less so it is in the solid phase. In fact, with 
JWF one gets solidification at a density much larger than the experimental 
value and the particles turn out to be very localized.^ A wave function that 
greatly improves over the JWF is the so called shadow wave function. In 
fact, SWF presently gives the best variational representation of the ground 
state properties of 4 He both in the liquid and in the solid phased 

In a SWF the atoms are correlated not only by standard direct cor- 
relations among particles, but also indirectly via the coupling to a set of 
subsidiary (shadow) variables.®! Integration over the shadows variables in- 
troduces effective interparticle correlations between pair of atoms but also 
between triplets and, in principle, to all higher orders in an implicit way; this 
is done so efficiently that it is possible to treat the liquid and the solid phase 
with the same functional formPSESI without the need to introduce a priori 
equilibrium positions to localize the atoms around lattice positions. Even 
for SWF in fact, solidification emerges as a result of a spontaneously broken 
translational symmetry due to the increased correlations as the density of 
the system increases. Moreover local disorder processes due to zero point 
motion are in principle allowed. This makes the use of such a wave function 
specially useful in the present context, where our purpose is to investigate 
the ground state nature of a quantum crystal. 

The general form of a SWF is given by 

= -^=MR) I dS G(R, S)MS), (3) 
v Qn J 

where R = (ri,...,rjy) and S = (si, . . . , sjv) are, respectively, the real 
and the shadow coordinates of the iV atoms and Qn is the normalization 
constant. As usual with SWF,® 1 4> r (R) is a Jastrow function which gives the 
direct correlations between particles, and we assume a McMillan form for 
the pseudopotentialP^ 



MR) = U e H?W • W 

i<j 

4> S (S) is another Jastrow function giving the inter shadow correlations; as 
pseudopotential we use the rescaled and shifted interatomic potential 

MS) = l[e- sv ^ . (5) 

i<j 



Q(R,S) is a Gaussian kernel coupling each shadow variable to the corre- 
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sponding real one: 

N 

S(R,S) = JJ e - c l r< - Si l 2 . (6) 

i=\ 

When SWF is employed as trial wave function in a variational descrip- 
tion of 4 He, the parameters b, m, 5, 7 and C have to be determined by 
minimizing the expectation value of the Hamiltonian 

h 2 N 

i=l i<j 

where V is the helium-helium interatomic potential. For instance when an 
hep commensurate crystal with N = 180 4 He atoms at p = 0.0293A -3 , 
slightly above the melting density, is considered and the standard Aziz 
potentiaP^is chosen as V, the minimization of gives the values b = 2.76A, 
m = 5, § = 0.1 IK" 1 , 7 = 0.875 and C = 0.8725A" 2 as optimal values for 
the variational parameters.^ Also SWF belong to the class of wave func- 
tions that display BEC even in the solid phase and a recent investigation 
computed the condensate fraction to be 5 x 10~ 6 at the melting density.^ 
When performing ground state calculations with a SWF by the Monte 
Carlo methods, due to the integration over shadow variables, two indepen- 
dent kinds of shadow particles must be taken into account in performing 
averages, as it can be easily seen, for example, by writing the expectation 
value of a diagonal observable 0{R) 

(O) = J dRdSdS'p(R,S,S')0(R) (8) 

where the probability density p(R, S, S') is given by 

p(r, s, s') = -^MS'MR, s')MR) 2 ®(R, s)MS) (9) 

Qn = J dRdSdS'cj) s (S')@(R, S')<fi r (R) 2 0(R, S)<j> s (S) . (10) 

A quantum-classical correspondence exists for SWF as well.™ In fact, 
there is a formal analogy between Eq. @ and the probability distribution 
density in configurational space of N classical triatomic molecules with flex- 
ible bonds at an inverse temperature (3 interacting via the potential 



U S (R,S,S') = ^ 



N 



i<j x ^ ' i=l 



+8Y,(V{ 1 s ij ) + V( ri s , ij )) 

i<j 
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Each molecule is then composed by one real (central) particle and two shad- 
ows. In this quantum-classical analogy no interaction other than the in- 
tramolecular forces corresponding to the real-shadow coupling and the inter- 
molecular forces between particle of the same species (provided by the real- 
real, shadow-shadow and shadow'-shadow' pseudopotentials) are present. 
The statistical sampling of ^(i?)) 2 then, maps the quantum system of N 
particles in an equivalent classical system of N flexible triatomic molecules 
with special interactions and the logarithm of the normalization constant 
(|10p is proportional to the excess free energy of such classical molecular 
system. Furthermore the Chester's argument can be extended also to this 
molecular solid to prove that also SWF gives rise to a finite equilibrium 
concentration of vacancies. 



3. EQUILIBRIUM CONCENTRATION OF VACANCIES 

As already pointed out, Chester's argumenlP^ ensures that in the quan- 
tum solids described by JWF and SWF (i.e. the quantum solids whose 
ground state wave functions are JWF and SWF) have a finite equilibrium 
concentration of vacancies, whose amount obviously depends on the chosen 
pseudopotential parameters. The problem of the estimate such equilibrium 
concentration of vacancies x v in the quantum solids at T = K is then 
formally equivalent to the computation of x v in their equivalent classical 
solids. The presence of a finite concentration of vacancies at equilibrium 
conditions means that the partition function Qtf of a macroscopic system of 
M particles has contributions from different pockets in configurational space 
corresponding not only to the commensurate state but also to states with a 
different number of vacancies; and the pockets with vacancies give the main 
contribution.^ This translates in the quantum case to say that the wave 
function ip{R) of a bulk system is describing at the same time both states 
with and without vacancies, and that the overwhelming contribution to the 
normalization constant Q_\f of the macroscopic system derives from pockets 
corresponding to a finite concentration of vacancies.^ 

Only about 10 years ago, x v has been computed for solid 4 He by exploit- 
ing this quantum-classical isomorphism^! anc l choosing as wave function a 
JWF with a McMillan pseudopotentiaP^. At melting density x v turned out 
to be x v = 6.36 x 10 -6 . This is interesting but one should keep in mind that 
JWF gives a very poor description of solid 4 HeJ 28 ^ 29 ^ moreover the pseudopo- 
tential parameters were chosen to stabilize the solid and not to minimize 
the expectation value of the quantum Hamiltonian at the chosen density. In 
addition, a quasi-harmonic approximation was used to estimate the free en- 
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ergy cost of a vacancy but the accuracy of this approximation is not known. 
Here we improve that estimate by employing direct simulations method to 
compute the involved free energies without any harmonic approximation. 
More important, with these methods we compute x v also for the SWF, a 
wave function that gives a very accurate description of solid 4 HeP^ 



3.1. Classical solid 

Pronk and FrenkeP^ outlined a full grand canonical route to estimate 
the equilibrium concentration of vacancies in the thermodynamic limit for a 
classical solid. The grand partition function Z of the crystal accounts for the 
fluctuations both of the number of particles N and of the number of lattice 
sites M. The canonical partition function of the crystal with n = M — N 
vacancies and M lattice sites is factorized, under the assumption of non- 
interacting defects, as 

Q "-^ T ^^w^ Q(n ' ivx > (12 » 

where Q^{V,T) is now the partition function of the crystal with n fixed 
vacancies at given lattice positions and the binomial factor accounts for all 
the possible configurations of such vacancies. This assumption allows to 
recognize in Z a binomial expansion over n that can be replaced with its 
sum value. The sum over the lattice site number M is approximated with a 
standard maximum term argument and then the obtained concentration of 
vacancies at the thermodynamic equilibrium is given by^l 

Xv = e -P(p-h) (13) 



where /i is the chemical potential of the defect-free crystal and —fx is the 
change in the free energy of the crystal due to the creation of a vacancy at 
a specific lattice position. 

At a first sight, the formula (jl3[) for x v seems to be obtained under the 
assumption of localized defect, and then should not be used in the case of 
mobile vacancies. This is not true, since the canonical partition function 
Qhi-n can be always factorized over the configurations of the defects even if 
mobile.^ In fact, the partition function Qm-u is made of static integrations 
over those particle configurations which are compatible with the presence 
of n vacancies. The corresponding volume in the configurational space can 
be split in sub-volumes each corresponding to a particular configuration of 
the n defects, i.e. Q^ n \ independently on how much time the defects spend 
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in any configuration. When vacancies are not interacting all these 



are 



equal and Eq. (|12p follows. 

Exploiting thermodynamic relations, the exponent in (|13p can be writ- 
ten as —/3(fi — fi) = (M— l)/3(/o — fd) — pP/p; then, in order to compute x v , 
we need the free energy per particle of the perfect, fo, and of the defected (i.e. 
with one vacancy), fd, crystal and of the pressure P. The pressure P for a 
classical solid is quite straightforwardly obtained with the virial method 3 - 3 -^ 
or from volume perturbations.^ Convergence of both the methods to the 
same value would provide a check on the estimated pressure. The computa- 
tion of the free energy (both /o an d fd) is less immediate: free energy cannot 
be directly obtained by a single Monte Carlo simulation.^ We recall that 
the free energy of the equivalent classical solid is related to the logarithm of 
the normalization constant Qn for the quantum crystal which, also, is never 
explicitly computed. 

Let us consider a classical system with interparticle potential energy 
U(R) where R denotes the set of coordinates. The thermodynamic integra- 
tion method gives a way for reconstructing free energy differences by inte- 
gration over a reversible path in the phase spaceP^ In a simulation we are 
not limited to use a thermodynamic path, rather all the parameters in the 
potential energy function can be used as thermodynamic variables. Then, 
in order to obtain the absolute free energy F, we must have at one end of 
the path a system whose free energy is known (reference system) , and to the 
other end we have the system of interest. For solids a popular method is the 
one of Frenkel and Lad d 36 * 3 \ The basic idea is to reversibly transform the 
solid under consideration into an Einstein crystal. To this end, the atoms 
are harmonically coupled to their lattice sites. Let U-Em(R) be the Einstein 
crystal chosen as reference system and let us consider the potential 



As A varies from to 1 the system described by IA\ goes from the system of 



U X (R) = (1- \)U{R) + XU Em (R) 



(14) 




and then, for the free energy 



F = F Ein + 



Jo 




dUx(R) 




where (• • • ) a denotes the canonical average for a system with potential energy 
function U\(R). It should be noted that the thermodynamic integration 
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method (|16p is intrinsically static, i.e. the derivative of the free energy is 
obtained in a sequence of equilibrium simulations each for different values of 
the coupling parameter A. The potential energy for the reference Einstein 
crystal is given by 

N 

^Ein(fi) = ^a i (r i -r°) 2 (17) 
i=i 

where r® are the lattice positions; its free energy is 




where, (3 is the inverse temperature. In the case of the defected crystal, 
a particle is removed keeping the simulation box volume and the lattice 
parameter unchanged. The spring constants can be adjusted to optimize 
the accuracy of the numerical integration in (|16p and are usually chosen to 
make the mean-squared displacement equal for A = and A = lP2 



3.2. Quantum solid 

In the quantum system U(R) corresponds to the so called pseudopo- 
tentials, for instance U(R) of the equivalent classical solid corresponds to 
log | ?/>j(i?) | 2 in the case of a JWF. The method outlined in the previous sub- 
section for computing x v might be used also for a quantum crystal described 
by a JWF or a SWF. However some modifications have to be introduced; due 
to the peculiar characteristics of the considered classical solids whose ther- 
mal equilibrium mimics the zero-point motion of quantum crystals, particles 
(or molecules) can be mobile. This is particularly true when a vacancy is 
present. In such situation the harmonic coupling to lattice positions would 
produce unbounded fluctuations because of the sudden suppression of ex- 
change processes, which are allowed by the wave function, when A ^ 0. A 
trick to bypass this problem is to break the thermodynamic path in two 
steps. In the first the equivalent classical crystal is changed into another 
crystal in which, by acting on the variational parameters, the pseudopo- 
tential are made so tight that the zero point motion is strongly reduced 
and exchange processes are greatly suppressed. In the second step this new 
crystal is transformed into the reference Einstein crystal constructed as pre- 
scribed by the Frenkel-Ladd methodP^ We have already underlined that in 
thermodynamic integration all the parameters in the potential function can 
be chosen as parameter to define the thermodynamic pathP^ 
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By construction the Frenkel-Ladd method gives access to Q("); in fact, 
the vacancy localization is explicitly introduced via the coupling to the Ein- 
stein crystal where it has a fixed position. This is true also for our two step 
thermodynamic path except than for the very beginning of the first step, 
where the vacancy is able to move. Strictly speaking this would introduce 
a systematic error, but in practice this error is smaller than the statistical 
one when only one vacancy is present. This has been verified by considering 
statistically independent estimates of (3fd which turn out to be compatible 
within the statistical error. In fact, thanks to the periodic boundary condi- 
tions the most of the sampled configurations (i.e. the ones where the vacancy 
occupies a lattice site) can be rigidly translated to keep the vacancy in the 
same original position without affecting the computed averages. 

We have checked our numerical code for the estimate of the free energy in 
classical systems with the Frenkel-Ladd method by reproducing the results in 
Ref. 1371 for a system of soft spheres, both with a single direct thermodynamic 
path from the interest crystal to the Einstein crystal, and with a two steps 
thermodynamic path as described above. It is known that the Frenkel- 
Ladd method is sensitive to the size scaling.^ We have considered then 
different boxes of increasing size (different number of lattice sites M) and 
we have found, for all the considered crystals the correct linear dependence 
of /3/+log M/M on l/MP^Since both (3fo and f3fd displays the same scaling 
law, the value of j3{p — /i), and then of x v , shows no appreciable dependence 
on M. 



4. RESULTS 

4.1. Equilibrium concentration of zero point vacancies for 

Jastrow wave function 

In a first computation we have considered the same wave function stud- 
ied in Ref. 1251 i.e. we have considered a JWF ([IJ) with a McMillan^ pseu- 
dopotential to describe a fee crystal in a box with M = 256 at p = 0.029A" 3 . 
The pseudopotential parameters are b = 4. 238 A and m = 6. As already no- 
ticed this b value does not minimize the expectation value of the Hamiltonian, 
but it is chosen so that the solid phase is stable at the considered density. 
To avoid any vacancy mobility problem when computing fd, the thermody- 
namic path has been split into two steps. The first step runs from the JWF 
crystal (i.e. the classical crystal equivalent to \tpj(R)\ 2 ) to a JWF' crystal 
obtained by increasing by a factor 1.4 the parameter b; this reduce the Lin- 
demann ratio from tl = 0.15 to tl = 0.05, which is small enough to ensure 
vacancy localization. The second step links this JWF' crystal to an Einstein 
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Fig. 1. Obtained values for (dj3U\(R)/d\)\ in a JWF fee crystal with M = 
256 at p = 0.029A~ 3 . Error bars are smaller than the used symbols, (a) first 
step: U X (R) = b x Uj{R) with Uj(R) as defined in Eq. ©, b x = (1 - X)b + 
X(b f /b) m and b f = 6.0A. (b) second step: U X (R) = (1- X)Uj'(R) + XU E m{R) 
where the harmonic constants of the Einstein crystal are chosen such that 
3/2/3a = 0.0318A 2 in the perfect crystal and 0.0355A 2 in the defected one. 



crystal of the type described by Eq. (|17|) with the same spring constant for 
all the particles. In order to ensure a good convergence of the results, we 
have considered at least 26 values of the coupling parameter for each step 
and for each value of the coupling parameter as many as 2.5 x 10 5 sweeps 
are produced at the equilibrium. The obtained results for (df3U\(R)/d\) \ 
are reported in Fig. [1] both for the first and the second thermodynamic path 
steps. The final result is — — /i) = —13.476 ± 0.057 which corresponds 
to an equilibrium concentration of vacancies x v = (1.4 ± 0.1) x 10~ 6 . Our 
result is about 5 times smaller than the previous estimate x v = 6.36 x 10~ 6 
of Ref. [25j This discrepancy is presumably due to the fact that in the older 
estimate, the free energy cost of a vacancy was computed relying on an 
harmonic approximation. Notice that the value of x v has an exponential 
dependence on —/?(// — /i), and it is easy to see that a value of — f3(/i — f\) 
11.2% higher than the one obtained in our computation would give the same 
x„ value as Ref. I 
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x v for JWF turns out to be very small. This is not too surprising 
because we know that a JWF can describe a solid but the particles turn out 
to be much more localized than in 4 He. 



4.2. Equilibrium concentration of zero point vacancies for 

shadow wave function 

The outlined method for the evaluation of the logarithm of the nor- 
malization constant Qn can be extended also to the molecular classical solid 
corresponding to the SWF. A convenient way to build the Einstein crystal for 
such SWF crystal is to couple each particle of a "molecule" (in the present 
case two shadows and one real) to the lattice position with an harmonic 
spring.— Since the choice of the pseudopotential parameters in a variational 
theory is related to a minimization principle, we choose as parameters the 
optimal ones to describe solid 4 He at p = 0.0293A -3 , a density which is 
slightly above melting in order to reasonably prevent instabilities. 

We have considered a simulation box which houses a hep crystal with 
M = 180 lattice positions and we have considered a two steps thermody- 
namic path. The first step starts from an optimized SWF (i.e. a crystal wave 
function where the variational parameters are set to minimize the expecta- 
tion value of the quantum Hamiltonian--) and ends on a SWF' crystal, ob- 
tained by increasing the real-shadow correlations (parameter C) by a factor 
6, and the shadow-shadow correlations (parameter 5) by a factor 4.5 in (jlll) . 
In this way the Lindemann ratio of the SWF crystal reduces to tl = 0.11 
from the starting tl = 0.27. The second step links this new SWF' crystal to 
the reference one: the Einstein crystal. Even in this case we have considered 
at least 26 values of the coupling parameter for each thermodynamic inte- 
gration step. Our results for (d(3U\(R)/d\) \ are reported in Fig. [2] for both 
the thermodynamic path steps. We have found — /3(/i — /i) = —6.56 ± 0.07; 
and then the estimated equilibrium concentration of vacancies for a SWF at 
the melting density is x v = (1.4 ± 0.1) x 10 -3 . 

We have checked our result by following different thermodynamic paths. 
We have considered different SWF' crystals, provided that the variational 
parameter were chosen in such a way to reduce the Lindemann ratio to values 
lower or comparable to those of standard classical solids. We have taken 
into account also two different kind of reference Einstein crystals: one with 
the same spring constant for all the particles and one with different spring 
constants for shadow and real particles. The same x v value is recovered in 
all the considered cases within error bars. 

The obtained x v for SWF is not so far from the actual experimental 
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Fig. 2. Obtained values for (d/3U\(R)/d\)\ in a SWF hep crystal with 
M = 180 at p = 0.0293A~ 3 . Error bars are smaller than the used symbols, 
(a) first step: 0U x {R) = (6/r^) m + C A £ti (|r 4 - s,| 2 + |r< - s',) 2 ) + 

Y.iKj (Vfrsij) + Ffrsy)) with Ca = (1 - A)C + AC/, C/ = 5.0A 2 , rJ A = 
(1 - A)r5 + X6 f and <5/ = 0.5K- 1 . (b) second step: ^ A (,R) = (1 - X)U S '(R) + 
A^Ein(-R) where the harmonic constants of the Einstein crystal are chosen 
such that 3/2/3a = 0.0766A 2 both in the perfect and in the defected crystal. 



bound to the zero-point vacancies concentration in 4 He of 0.4%.-^ Then 
improvements in experimental investigations could answer the question if 
SWF describes well also the equilibrium concentration of vacancies in the 
ground state. Moreover, by using the BEC transition temperature Tbec for 
an ideal Bose gas with mass equal to the vacancy effective mass as obtained 
with SWF^Sl (m v = 0.35mif e ) as an estimate of the supersolid transition 
temperature, we obtain T$s — ll.Sxy = 141 mK, which is about 2.3 times 
larger than the experimental estimated transition temperature T = 60 mK. 7 
One possible origin of this discrepancy is the fact that x v is estimated via 
a variational technique. As we shall discuss later, the outlined route to 
compute x v can be in principle extended to exact techniques. Another origin 
is the presence of a significant vacancy-vacancy interaction. 
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Fig. 3. Vacancy- vacancy correlation function g vv (r) computed in an hep crys- 
tal at the melting density containing three vacancies and M = 540 lattice 
positions. g V v( r )/ grand is the vacancy- vacancy correlation function normal- 
ized with the correlation function of an ideal gas of 3 particles on the same 
hep lattice, while in the inset g vv (r) is reported normalized to unity. 



4.3. Vacancy vacancy correlations 

Eq. (113j) for x v is valid under the assumption of non interacting defects. 
A key aspect to be investigated is whether in the quantum solid described by 
SWF vacancies interact and then if a dilute gas of vacancies is stable. There- 
fore we have studied systems with two and three vacancies. Here we report 
the results for a system of N = 537 atoms (plus twice as much shadows) at 
density p = 0.0293A" 3 in a box with periodic boundary condition such that 
an hep crystal with M = 540 lattice positions fits in. Starting from an ini- 
tial configuration corresponding to such hep lattice in which 3 particles have 
been removed, we find that the crystalline state is stable and the periodic- 
ity is such that the number of density maxima is equal to 540, i.e. 3 more 
of the 537 atoms. This means that we have an hep crystal with 3 mobile 
vacancies. We have monitored the 3 vacancies positions during the Monte 
Carlo sampling, and by histogramming their relative distances on at least 
4 x 10 6 equilibrium configuration we have constructed a vacancy-vacancy 
correlation function g vv (r). In the inset of Fig. [3l we plot the obtained 
<7v V (r) normalized to unity, which gives the probability, given one vacancy at 
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the origin, to find another one at distance r. In the main of Fig. [3] we report 
g vv (r) normalized with the correlation function of an ideal gas of 3 particles 
on the same hep lattice. In this case we have a vacancy pair distribution 
function. The determination of the vector positions of the vacancies in a 
crystalline configurations is far from being trivial due to the large zero point 
motion, to the vacancy being mobile and because the center of mass is not 
fixed. Our algorithm proceeded as follows: given a configuration explored 
by the Metropolis sampling, we find the crystal lattice which best fits the 
positions of the real particles; this is achieved through the construction of 
a "Gaussian local density" model which is optimized adjusting the position 
of the center of mass of the lattice so that the "molecules" (a real and two 
shadow variables) fall in high density regions. Then, the position of a va- 
cancy is defined as that site which does not host any atom in its Wigner-Seiz 
cell and whose nearest neighbors are not double occupied!^ 

From the plot of g vv (r) (Figj3]), we see that the vacancies explore the 
whole available distance range. We may so conclude that no sign of bound 
state is detected, at least one the scale of the simulation box, and it seems 
quite improbable a larger range bound state. Similar results is found for two 
vacancies. The observed large value for g vv (r) at nearest neighbor distance 
is an indication of some short range attraction. But this interaction seems 
not large enough to give a bound state at least for up to 3 vacancies, as 
evident in the tail of <7w( r )/#rand which approaches the unity, i.e. vacancies 
behave almost like ideal particles for distances beyond about 10 A. We then 
conclude that the vacancy gas should be stable in bulk SWF solid; and, for 
very low concentrations, as the expected equilibrium one, the approximation 
of vacancies as non interacting defects seems to be reasonable. However, the 
effect of an attractive interaction between vacancies, if not strong enough 
to lead to phase separation, would affect the equilibrium concentration of 
vacancies. For example, the concentration of divacancies is given by 



X2v - Q-PAf 



where q is the coordination number of the considered lattice {q = 12 for hep 
and fee) and A/ is the gain in the free energy in bringing two single vacancies 
together to form a pair!^ The relation (|19j) is the one expected from the 
application of the mass action law to the vacancy association reaction!^ The 
overall effect is that the total concentration of vacancies should somewhat 
increase with respect to the value computed in the previous Section. 



(19) 
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5. CONCLUSIONS 



We have presented the first quantitative estimate of the concentration 
of vacancies for a quantum solid described by the SWF which has been vari- 
ationally optimized to describe solid 4 He. We have computed x v exploiting 
the formal analogy between the ground state of bosons and the Boltzmann 
weight for classical particles and using advanced simulation methods. At 
melting we find x v = (1.4 ± 0.1) x 10~ 3 . This is a sizable concentration not 
too far from the present experiment upper bound for zero point vacancies in 
solid 4 He. x v for SWF is much larger than the value x v = (1.4 ± 0.1) x 10~ 6 
given by a JWF. This, presumably, is due to the much stronger localization 
of particles described by a JWF compared to SWF. The study of two and 
three vacancies in the hep crystal within the SWF shows the presence of a 
significant short range attraction but such multiple vacancies do not form a 
bound state. 

SWF has been proved to be so successful in describing many properties 
of solid 4 He that it actually represents the best available variational descrip- 
tion of 4 He. As all the variational descriptions it suffers of some limitations, 
for example SWF gives a finite condensate fraction also for a commensurate 
crystal which turns out to be zero when computed with exact techniques 
such as PIMC at finite temperature™!] and SPIGS at zero temperature.- 
It has been suggested that this defect of SWF might be due to a qualitative 
lack in the functional form of the wave function, which, probably, misses 
the description of zero point motion of excitations other than longitudi- 
nal phonons, such as the transverse onesP^For this reason some prudence is 
needed with respect to the relevance of the obtained concentration of ground 
state vacancies for SWF to the real ground state of solid 4 He. For example, 
in Ref. [TO] it was argued that if the true ground state has a finite concen- 
tration of vacancies, the energy per particle computed in a box with the 
correct x v should be lower than the one computed for the commensurate 
crystal. When the Aziz potential^ is used as interaction potential for solid 
4 He in the quantum Hamiltonian ([7]), this property is unfulfilled by SWF 
at their equilibrium concentration of vacancies. This will be true only if 
the energy per particle is computed with the potential that has the SWF as 
exact ground state in ([7]). Unfortunately such a potential is unknown; while 
it is known the one that has the JWF as exact ground stated It is easy to 
check, as pointed out on a general discussion in Ref. [241 that this is not a 
simple pair potential but it has a very strong 3-body term, which is about 
1.3 times the 2-body one in the considered case. 

The outlined route for the estimate of the equilibrium concentration of 
vacancies can, in principle, be extended also to the exact SPIGSp^ method 
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which recovers more and more accurate approximations of the true ground 
state wave function by successive short time projections in imaginary time 
starting from a SWF. The classical solid equivalent to SPIGS turns out 
to be a collection of open polymers whose length (number of constituent 
monomers) is fixed by the number of considered projection steps. The 
Chester's argument applies also to this polymeric crystal which is expected 
to have an equilibrium concentration of vacancies. Further investigation 
with exact technique is not trivial but advisable. 

This work was supported by the INFM Parallel Computing Initiative, 
by Supercomputing facilities of CILEA and by Mathematics Department "F. 
Enriques" of the Universita degli Studi di Milano. 
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